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The distributions of work for strongly non-equilibrium processes are studied using a very general form of a 
large-deviation approach, which allows one to study distributions of almost arbitrary quantities of interest for 
equilibrium, non-equilibrium stationary and even non-stationary processes. The method is applied to varying 
quickly the external field in a wide range B = 3 o for critical (T = 2.269) two-dimensional Ising system of 
size L x L = 128 x 128. To obtain free energy differences from the work distributions, they must be studied in 
ranges where the probabilities are as small as 1CF 240 , which is not possible using direct simulation approaches. 
By comparison with the exact free energies, which are available for this model for the zero-field case, one sees 
that the present approach allows one to obtain the free energy with a very high relative precision of 1(T 4 . This 
works well also for non-zero field, i.e., for a case where standard umbrella-sampling methods seem to be not so 
efficient to calculate free energies. Furthermore, for the present case it is verified that the resulting distributions 
of work for forward and backward process fulfill Crooks theorem with high precision. Finally, the free energy 
for the Ising magnet as a function of the field strength is obtained. 
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Studying non-equilibrium work processes [Q] 01 has be- 
come a useful tool to extract information about physical sys- 
tems. Particular useful are the Jarzynski relation y3J and the 
related Crooks theorem fl, which allow one to extract equi- 
librium information from non-equilibrium systems. Here, a 
system in contact with a heat bath at temperature T is prepared 
initially in equilibrium, under the influence of some external 
parameter B = B\. Next it is driven out of equilibrium via 
quickly [5J] varying B = B\ — >• B2, into another state, while 
performing some work W. The Jarzynski relation relates the 
free energy difference AF between the equilibrium states at 
B = B\ and B = B2 to the work W performed during the 
non-equilibrium process B = B\ — >• B^: 



-AF/T 



(1) 



where (...) denotes the combined average over the initial 
equilibrium distribution and all possible process paths. In a 
similar way, the theorem of Crooks relates the distributions of 
work P(W) for the forward process and of the negated work 
P ICV (—W) for the reverse process B = B2 — > B\ (where one 
starts in equilibrium at B = B2) to the same AF via: 



P(W) 

Pr CV (-W) 



(W-AF)/T 



(2) 



Hence, P(W) and P tev (-W) should intersect at W* = AF. 
Unfortunately, the average of ([TJ and the region where the two 
distributions intersect are often dominated by exponentially 
small probabilities, making finite-sampling errors particular 
strong Q6Q. Thus, the author of this work is only aware of ap- 
plications which exhibit either a small number of degrees of 
freedom, or where the initial and final states B = B\ , B2 are 
very similar to each other. E.g., for the Ising model in an ex- 
ternal field B, work distributions have been obtained iM, 



for rather large systems exhibiting N = 128 2 spins, but suc- 
cessfully only in the paramagnetic and in the ferromagnetic 
phases, where the work distribution could be approximated 
very well by a Gaussian. At the critical point, it was only 
possible to sample the work distribution reliably for very slow 
and small changes of the field, making the direct application 
of the approach not successful. 

As it is shown in this work, the work distributions can 
be obtained very reliably via large-deviation or importance- 
sampling techniques ]9[], which are able to address large- 
deviation regions lioll of interest using bias functions. Such 
techniques have been previously applied numerically to study 



work distributions for small systems 111 ll 11211 . If one targets 



not obtaining the full work distribution but is interested just 
in free ener gy d ifferences, it was suggested from results of 
simulations 1 1130 and analytical studies 111411 of small model 
systems, that applying work theorems cannot compete with 
direct umbrella-sampling techniques which explicitly obtain 
the distribution of the energy over large ranges of the support. 
Thus, these results were rather discouraging. Nevertheless, in 
the present work not only work distributions are obtained over 
a large range of the support, but it is also shown for a sam- 
ple strongly non-equilibrium process in a large system with 
a non-Gaussian work distribution that a large-deviation non- 
equilibrium work-sampling approach turns out to give very 
accurate results. 

The algorithm presented in this work is a very general 
"black-box" type approach which renders it applicable to 
study the distribution of almost any quantity of interest 
for equilibrium, non-equilibrium stationary and even non- 
stationary processes. The algorithm is here applied to work 
distributions of the Ising model in a non-zero field. In previ- 
ous work the free energy could be obtained using umbrella- 
sampling approaches for only rather small systems 111 51 . This 
is in contrast to the zero-field case, where indeed umbrella 
sampling is most efficient. Here, the work distributions for 



the non-zero-field case are directly obtained down to proba- 
bilities as small as 10~ 240 such that the Jarzynski relation ([T) 
and Crooks theorem © can be directly evaluated. For this 
purpose, the explicit biased sampling not only over the paths 
but also over the initial equilibrium distribution is included 
and simulations are performed for a large range of freely 
adjustable bias weights, allowing do obtain the distributions 
over hundreds of decades in the probability. The simulations 
are performed for large systems and strongly non-equilibrium 
paths. 

The ferromagnetic Ising model in a field B > is studied 
here, given by a set of TV Ising spins s; = ± 1 and described 
by the Hamiltonian H = —JJ2(ij) s i s j ~~ ^J2i s i- The 
first sum runs over all bonds connecting neighboring sites of 
a square lattice of size L, i.e., N = L x L. Periodic boundary 
conditions are applied in both directions. The system is cou- 
pled to a heat bath at temperature T = 2.269, about the critical 
temperature for the ferromagnet-paramagnet phase transition. 

Next, the numerical approaches are described. Processes 
were considered, where the system is started in equilibrium 
and within niter steps the field is changed. For the forward 
process it was increased from B = to B = -B max , i-e-, 
in each step by AB = B max /nit e x- Thus, during each 
step I the work Wi = — ^ SiAB was performed, the to- 
tal work is W = Yli Wi- After each field increment, one 
sweep of a Monte Carlo (MC) simulation llql with single- 
spin flip Metropolis dynamics was performed. Hence, in each 
sweep, N times a spin was randomly chosen and a spin flip, 
exhibiting an energy change AH, was accepted with prob- 
ability min{l, exp(— AH/T)}. The initial equilibrium con- 
figuration was obtained by starting in a random configuration 
and performing 1000 steps of the Wolff cluster algorithm II 1711 . 
which should ensure equilibration since the auto-correlation 
time for this algorithm is of the order of r 10 at T c 1 1811 . 
Also for B max = 3 the reverse process was considered where 
the system was started in equilibrium at B = S max and the 
process B = B max — > is performed in an analogous way as 
for the forward process. Since the equilibrium configuration 
for this case is almost fully magnetized at this large value of 
B = B max (typically 0.002A spins are flipped), the initial 
configurations were obtained by starting with all spins up and 
performing one sweep of the single-spin flip dynamics prior 
to the B = _/3 max — » process. 

For the case A = 128 2 and 10 6 independent simulations, 
the histograms P(W) of work for the forward and P lev (— W) 
for the reverse process are shown in Fig. [T] According to 
the theorem of Crooks, the two histograms should intersect 
at W* = AF. This appears to be somewhere between 
W = -50000 and W = -45000, which on the first sight 
is only a small interval compared the support of the distri- 
bution P(W) visible in Fig.Q] Nevertheless, the probabili- 
ties become so small, that it is impossible to see the intersec- 
tion using any feasible number of standard simulations runs. 
Actually, as it is shown below, the crossing appears where 
P(W) = P iev (—W) ~ 10~ 57 . Hence, numerical large devi- 
ations techniques have to be used, to address this region. 
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FIG. 1: Simple-sampling distribution P(W) of work for an Ising 
system with N = 128 2 spins for the forward 5 = 0—5-3 and the 
mirrored distribution P rC v( — W) for the reverse process B = 3 — > 0. 



The algorithm presented here is different compared to well- 
known algorithms as, e.g., the "cloning" approach IU9H21I1 . 
and consists of a second MC-simulation level. Each config- 
uration is represented by a vector £ = (£1,^2, ■■ • ,£,m) of 
suitable length M (see below). The basic idea is that the 
entries of £ are random variables uniformly distributed in 
[0, 1], which are used to feed the random process under in- 
vestigation. Hence, e.g., when performing the work distribu- 
tions, the random decisions are not based on numbers drawn 
from random number generators, but, in a defined manner, 
on the entries of £. Here, for the forward process, the first 
«Woiff(2A + 1) entries are used to feed nwoiff = 10 iter- 
ations of the Wolff algorithm, starting from a precomputed 
(1000 Wolff iterations) equilibrium configuration s^-°\ This 
allows to sample the equilibrium distribution prior to the work 
process. Each Wolff iteration consist of choosing one seed 
spin (consuming one entry of £) plus a cluster growth, where 
possibly for each of the 2N bonds it has to be decided ran- 
domly whether it is "activated" or not, for details of the Wolff 
algorithm see Ref. 0. Note that each bond is assigned 
a specific entry of £ (for each Wolff iteration, respectively), 
independent of whether the bond is tried to be activated or 
not. Next, the niter work sweeps are performed, consisting 
of njtcr — 1 single-spin-flip MC sweeps (the last sweep af- 
ter the final field increment can be omitted), where for each 
sweep 2N entries of £ are consumed, one for randomly se- 
lecting a spin, and one for the Metropolis criterion (also if 
AH < 0). Hence, for one full process M = ny? Q ig(2N + 
1) + ("iter — 1)2 A entries, corresponding to random num- 
bers, are used. For the reverse process, where no Wolff algo- 
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rithm is used, M — 2N + (m ter - l)2N = n itol -2N. Each 
process fed by £ results in a work W(£). Now, the second- 
level MC consist of changing a small number (here 400) 
of the entries of £, each drawn again uniformly from [0, 1], 
leading to a "trial configuration" £' with corresponding work 
W(£'). The new configuration is accepted with the Metropo- 
lis probability min{l, exp(-(W(£') - W(£))/T M c)}, other- 
wise £ is kept for the next second-level MC step. Thus, the 
observed distribution of work will exhibt a Boltzmann bias 
p T MC (W) ~ P{W)e- w/Tmc . Note that T MC is a freely ad- 
justable parameter, different from the temperature T of the 
Ising system, which allows to center the observed distribution 
in different regions. Therefore, by performing the simulation 
for different values of Tmc sucn that the resulting distribu- 
tions Pt mc (W) overlap, one can reconstruct the desired dis- 
tribution P(W) via reweighting and gluing the distributions 
together 112211 . Note that instead of using a Boltzmann bias for 
the observed work, one could also use an umbrella sampling 
on the £ vectors to obtain P(W) directly. However, this was 
tried during the present work extensively, but the Boltzmann 
bias was more efficient. 

Seen in an abstract manner, the approach is based on a 
configuration vector £, an evaluation function a sim- 

ple dynamics changing the configuration vector to create trial 
vectors £' and a Metropolis criterion (involving "tempera- 
ture" Tmc), which accepts £ based on H(£) and &(£'). 
All problem-specific details are included in the function H. 
Hence, if H was the energy of a spin system, then the second- 
level MC would be a standard MC simulation. But H can 
represent almost any process, like for the present application 
where it states the work W arising from an equilibrium sam- 
pling followed by a strongly non-equilibrium work process. 
Using this general view on the large-deviation simulation, it 
is clear that the approach can be applied basically to any equi- 
librium, non-equilibrium stationary and even non-stationary 
process which transforms a vector of a random numbers into 
a measurable result, independent of how involved this trans- 
formation is. Hence, the approach should be applicable to any 
random process which can be simulated on a computer. 

For the simulations, a number between 4 (£> max = 0.25) 
and 37 (T/ max = 3) of temperatures Tmc £ [1.5,200] were 
considered. The spacing ATmc between the MC tempera- 
tures ranged from 0.1 for low values of Tmc u P to 100 for 
larges values. For each value of Tmc. 10 6 MC trials were per- 
formed, taking about 6 hours on a core of standard 2.66 GHz 
Intel Westmere processor, i.e., just 37 x 6 = 222 core hours 
for the strongest field T max = 3. 

Concerning the analysis of the simulation leading to the full 
work distribution, the intersection region of the distributions 
of work for the forward and reverse processes are shown 
in Fig. |2] for B max = 3. The two distributions P(W) and 
P Iev (—W) intersect at W* re -47443. For comparison, 
also the exact free energy difference was obtained. For the 
zero-field case, the exact free energy To is known analytically 
1 2311 for finite-size systems. For the case B > 0, if B is large, 
the system is almost fully magnetized, except for a few single- 
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FIG. 2: Distribution of works for the forward and mirrored distribu- 
tion for the reverse process for B — o 3, L — 128. According to 
the theorem of Crooks, these distributions are expected to intersect 
at W* = AT signifying the free-energy difference AT. The ex- 
act free-energy difference is indicated by the vertical line. The inset 
shows a blow up of the intersection region. 



spin excitations. Hence, the free energy is given via Fb 



-T 



log (f 



(2+B)N/T 



where only few terms of the sum around the typical number 
k of excited spins (about 0.0027V for B = 3) contribute 
significantly. For L = 128, B = 3.0, this results in 
AT = Fb — To re -47433. Thus, the relative deviation of 
the estimated from the exact free energy difference is only 
(AT - W*)/AF = 0.0002. 

To verify whether the data fulfills Crooks theorem (f2), the 
histogram for the reverse process was rescaled accordingly, 
see Fig. [3] This is confirmed by the data with high precision. 
Note that testing Crooks theorem may also conveniently serve 
as a check that the second-level MC simulations are equili- 
brated. 

To obtain AT using the Jarzynski relation (|T), the integral 
(exp(-W/T)) = / exp(-W/T) P(W) dW has to be eval- 
uated, resulting in AT re —47438, which has a relative de- 
viation 0.0001 from the exact result. The integrand is shown 
in Fig. @] Note that only the region close to the peak around 
W = —48100 contributes significantly to the integral, which 
deviates much from the point where P(W) and P re v{—W) 
intersect. Nevertheless, one has to obtain the full distribution 
in a region ranging from its peak value at W = —42000 down 
to W = —48200 to get P(W) right. For the reverse process, 
the evaluation (see inset of Fig. [3]) results in — AT = 47450, 
which deviates by a factor of 0.0004 from the exact value. 
Note that one has to obtain T rcv (W) over an even broader sup- 
port, which is probably the reason for the somehow smaller 
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FIG. 3: Distribution of works for the forward and rescaled-mirrored 
distribution for the reverse process. The inset exhibits the integrand 
P TCV (W)exp(-W/T) used to obtain (exp(-W/T)) from the re- 
verse process B — 3 — > 0. 
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FIG. 4: Integrand P(W) exp(-W/T) used to obtain 
(exp(— W/T)}. Inset: Resulting free energy difference AF 
per spin N (minus Bmax) as a function of the final magnetic field 
Bmax for the forward process. 



accuracy compared to the forward process. 

The forward process was performed for different values of 
B max , see inset of Fig. [4] The amount by which AF/N is 
larger than B max describes the entropy loss due to the align- 
ment of the spins to the field. 

To summarize, a biased sampling approach is introduced, 



which is based only on a Markov-chain evolution of a vector 
of entries from the interval [0,1], seen as an input vector of 
random numbers to an arbitrary stochastic process, which can 
be treated as a black box within the approach. 

Here, high-precision determination of work distributions 
for Ising magnets in a field were performed, for large sys- 
tems and strongly non-equilibrium processes, hence for cases 
where traditional direct approaches for measuring work distri- 
butions completely fail. In the past only close-to equilibrium 
processes could be studied with small accuracy 0.0] ■ Still, 
the path sampling applied here is very general, since no details 
of the path construction have to be known to efficiently sam- 
ple the corresponding work distribution down to probabilities 
as small as 10 -240 . This contrasts the approach with previ- 
ous problem-specific yet quite successful techniques, like the 
"shooting approach" 112411 or "cloning" IU9U21I1 . 

The work distributions are used to extract free energy differ- 
ences for Ising systems in a field using the Jarzynski relation 
H and the theorem of Crooks fl. Note that for determin- 
ing the free energy of the Ising system without a field, very 
good other approaches exist. Using the convenient Wang- 



Landau umbrella sampling II25I1 . the free energy was deter- 
mined very accurately for systems of size TV = 256 2 . Based 
on measuring the large-deviation properties of the number 
of components for Fortuin-Kasteleyn clusters even systems 
of size TV = 1000 2 could be treated JH]. Recently it was 



claimed 111311 that in general umbrella-sampling approaches 
should be superior or at least equally efficient as applying 
large-deviation techniques to work distributions to measure 
free-energy differences. Nevertheless, for the present study 
of an Ising systems in a field, using umbrella-sampling only 
sizes of N = 42 2 could be studied so far |15 i|, about ten times 
smaller than the sizes addressed in the present study. Hence, 
for certain systems, e.g., for an Ising magnet in a field, al- 
most the full work distribution can be determined using the 
present approach. But even aiming only at determining free 
energy differences, the present very general approach might 
be superior to highly-evolved existing techniques. The ob- 
served failure of Refs. II13[ Il4ll might be due to the fact that 



there only very small systems could be studied. Also it could 
be due to the specific "shooting" algorithm 12411 which might 
not be most suitable for fully random processes. Finally, for 
some past studies also single IU2I1 or very many 111 ill specific 
bias functions where used, while here a Boltzmann reweight- 
ing for few selected weights was performed, which allows to 
address different regions of interest independently. 

Due to the black-box structure of the algorithm presented 
here, it allows to study equilibrium, non-equilibrium station- 
ary and even non-stationary systems. Hence, for future work, 
many applications of this algorithm can be anticipated. 
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